setwd("~/Google Drive/MDCurso/Datos")
datos<-read.csv("Ej1kmeans.csv",sep = ";",header=F)
head(datos)## V1 V2
## 1 -0.3508569 0.3348495
## 2 0.4312463 -0.2319262
## 3 -0.2342527 -0.1372323
## 4 -0.2942210 -0.4817277
## 5 -0.1512742 0.3637852
## 6 0.3247037 -0.1647549
grupos<-kmeans(datos,centers=2,iter.max=100) ## iter.max por defecto es 10
grupos## K-means clustering with 2 clusters of sizes 51, 49
##
## Cluster means:
## V1 V2
## 1 0.02169424 0.08865999
## 2 0.99128291 1.07898833
##
## Clustering vector:
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2
## [71] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##
## Within cluster sum of squares by cluster:
## [1] 7.489019 9.397754
## (between_SS / total_SS = 74.0 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
grupos$cluster## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2
## [71] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
grupos$centers## V1 V2
## 1 0.02169424 0.08865999
## 2 0.99128291 1.07898833
grupos$totss # Inercia Total## [1] 64.88887
grupos$withinss # Inercia Intra-clases por grupo (una para cada grupo)## [1] 7.489019 9.397754
grupos$tot.withinss # Inercia Intra-clases## [1] 16.88677
grupos$betweenss # inercia Inter-clases## [1] 48.0021
# Verificación del Teorema de Fisher
grupos$totss==grupos$tot.withinss+grupos$betweenss ## [1] TRUE
grupos$size # Tamaño de las clases## [1] 51 49
plot(datos,pch=19,xlab=expression(x[1]),ylab=expression(x[2]))
points(grupos$centers,pch=19,col="#FF9C5B",cex=2)
points(datos,col=grupos$cluster+1,pch=19)library(cluster)
library("factoextra")Loading required package: ggplot2
Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
datos.escalado <- data.frame(scale(datos), grupos$cluster)
clusplot(datos.escalado, grupos$cluster, color=TRUE, shade=TRUE,
labels=2, lines=0)barplot(grupos$centers[1,],col='#F5634A')barplot(grupos$centers[2,],col='#ED303C')rownames(grupos$centers)<-c("Cluster 1","Cluster 2")
barplot(t(grupos$centers),beside=TRUE,col=c("#F5634A","#F5634A","#ED303C","#ED303C"))## con iter.max=10 algunas veces encuentra soluciones suboptimales
grupos<-kmeans(datos,centers=4,iter.max=50)
grupos$cluster## [1] 1 4 2 2 1 4 1 2 4 1 1 1 4 4 2 4 1 2 2 4 4 4 4 4 2 4 2 2 4 2 4 2 2 2 2
## [36] 2 4 2 2 4 2 1 4 2 1 1 4 1 2 2 3 3 3 3 3 3 1 4 3 3 3 3 3 3 3 3 3 3 3 3
## [71] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
grupos$centers## V1 V2
## 1 -0.06821724 0.50099826
## 2 -0.19304414 -0.02509277
## 3 1.00982118 1.07626940
## 4 0.32002034 0.01295400
plot(datos,pch=19,xlab=expression(x[1]),ylab=expression(x[2]))
points(grupos$centers,pch=19,col="black",cex=2)
points(datos,col=grupos$cluster+1,pch=19)datos.escalado<- data.frame(scale(datos), grupos$cluster)
clusplot(datos.escalado, grupos$cluster, color=TRUE, shade=TRUE,
labels=2, lines=0)setwd("~/Google Drive/MDCurso/Datos")
datos <- read.csv("EjemploClientesCorregidaEdad.csv",header=TRUE, sep=";", dec=",", row.names=1)
str(datos)## 'data.frame': 37 obs. of 12 variables:
## $ Edad.10 : num 2.5 2.4 2.8 2.3 4.9 3.2 2.6 2.3 2.5 2.9 ...
## $ Antiguedad : int 1 0 7 0 6 4 0 4 4 0 ...
## $ Espacios.Parqueo : num 7.6 4.8 6.8 3.4 7 5.6 6.2 5.6 4.6 5.4 ...
## $ Velocidad.Cajas : num 7.6 9 8.4 7.8 3.2 7.8 8 6.8 8 6.4 ...
## $ Distribucion.Productos: num 7.8 7.2 7.6 9 1.2 6.8 6.6 6.2 3.8 8.8 ...
## $ Atencion.Empleados : num 9.7 10 8.7 10 10 10 9.3 9.7 10 9.7 ...
## $ Calidad.Instalaciones : num 5 2 2.7 1 4 3 3.3 4 1.7 6.7 ...
## $ Ubicacion : num 9 9.6 9.2 10 9 10 8.6 6.8 9.8 10 ...
## $ Limpieza : num 7.6 6.8 6.2 4.4 1.4 5 7.8 6.8 5 5.6 ...
## $ Variedad.Productos : num 5.6 8.4 9 4 4.8 4.2 6.4 7.4 4.4 6.2 ...
## $ Prestigio.Empresa : num 7 9.8 9.6 2.8 2.6 4.2 9.6 5.6 6 8.4 ...
## $ Calidad.Servicio : num 6.6 5.4 8.5 5.4 3.3 7.2 6.5 4.5 7.6 6.5 ...
grupos<-kmeans(datos,centers=3,iter.max=200)
grupos$cluster # Cluster al que pertenece cada fila## Ariana Guiselle Francisco Griselda Damaris Johana
## 3 3 1 3 2 1
## Bernal Freddy Estafania Laura Arnoldo Beatriz
## 3 1 1 3 1 1
## Rebeca Sofia Ingrid Rocio Karen Luis
## 3 1 1 1 1 1
## Pedro Lorena Elena Julian Natalie Shirley
## 3 1 1 2 3 3
## Andres Alejandro Grace Nuria Flor Roberto
## 1 1 3 3 3 1
## Victor Arturo Maritza Diana Juan Guillermo
## 3 1 1 3 1 3
## Silvia
## 3
grupos$centers # Centros de los clústeres## Edad.10 Antiguedad Espacios.Parqueo Velocidad.Cajas
## 1 3.057895 5.210526 6.4 7.926316
## 2 4.450000 5.500000 6.0 5.500000
## 3 2.475000 0.687500 5.8 8.162500
## Distribucion.Productos Atencion.Empleados Calidad.Instalaciones
## 1 6.389474 9.615789 3.405263
## 2 2.300000 9.850000 3.650000
## 3 7.237500 9.650000 3.318750
## Ubicacion Limpieza Variedad.Productos Prestigio.Empresa Calidad.Servicio
## 1 8.694737 6.663158 6.905263 6.894737 5.642105
## 2 9.100000 2.800000 5.400000 3.800000 2.400000
## 3 9.350000 7.200000 6.962500 7.700000 4.812500
colnames(datos) ## [1] "Edad.10" "Antiguedad"
## [3] "Espacios.Parqueo" "Velocidad.Cajas"
## [5] "Distribucion.Productos" "Atencion.Empleados"
## [7] "Calidad.Instalaciones" "Ubicacion"
## [9] "Limpieza" "Variedad.Productos"
## [11] "Prestigio.Empresa" "Calidad.Servicio"
color <- c("#FF6449", "#FEB035", "#FCE020", "#F7EC69", "#F1F8BE","#D5B9F6",
"#A2E3CD","#EDF380", "#D8FD9C", "#AEEC64", "#F598F8", "#9EFF37")
barplot(grupos$centers[1,],col=color,las=2,cex.names = 0.7, ylim = c(0,12))barplot(grupos$centers[2,],col=color,las=2,cex.names = 0.7, ylim = c(0,12))barplot(grupos$centers[3,],col=color,las=2,cex.names = 0.7, ylim = c(0,12))rownames(grupos$centers)<-c("Cluster 1","Cluster 2","Cluster 3")
barplot(t(grupos$centers),beside=TRUE,legend=colnames(datos),main = "Gráfico de Interpretación de Clases",col=color,ylim=c(0,30),cex.names = 0.7)# Gráfico Tipo Araña
library(fmsb)
centros<-grupos$centers
rownames(centros)<-c("Cluster 1","Cluster 2","Cluster 3")
centros<-as.data.frame(centros)
maximos<-apply(centros,2,max)
minimos<-apply(centros,2,min)
centros<-rbind(minimos,centros)
centros<-rbind(maximos,centros)
centros## Edad.10 Antiguedad Espacios.Parqueo Velocidad.Cajas
## 1 4.450000 5.500000 6.4 8.162500
## 11 2.475000 0.687500 5.8 5.500000
## Cluster 1 3.057895 5.210526 6.4 7.926316
## Cluster 2 4.450000 5.500000 6.0 5.500000
## Cluster 3 2.475000 0.687500 5.8 8.162500
## Distribucion.Productos Atencion.Empleados Calidad.Instalaciones
## 1 7.237500 9.850000 3.650000
## 11 2.300000 9.615789 3.318750
## Cluster 1 6.389474 9.615789 3.405263
## Cluster 2 2.300000 9.850000 3.650000
## Cluster 3 7.237500 9.650000 3.318750
## Ubicacion Limpieza Variedad.Productos Prestigio.Empresa
## 1 9.350000 7.200000 6.962500 7.700000
## 11 8.694737 2.800000 5.400000 3.800000
## Cluster 1 8.694737 6.663158 6.905263 6.894737
## Cluster 2 9.100000 2.800000 5.400000 3.800000
## Cluster 3 9.350000 7.200000 6.962500 7.700000
## Calidad.Servicio
## 1 5.642105
## 11 2.400000
## Cluster 1 5.642105
## Cluster 2 2.400000
## Cluster 3 4.812500
color <- c("#FCEBB6","#78C0A8","#5E412F")
radarchart(as.data.frame(centros),maxmin=TRUE,axistype=4,axislabcol="slategray4",
centerzero=FALSE,seg=8, cglcol="gray67",
pcol=color,plty=1,plwd=5,title="Comparación de clústeres")
legenda <-legend(1.5,1, legend=c("Cluster 1","Cluster 2","Cluster 3"),
seg.len=-1.4,title="Clústeres",pch=21,bty="n" ,lwd=3, y.intersp=1,
horiz=FALSE,col=color)# En grupos$cluster está a qué cluster pertenece cada fila de la tabla de datos
NDatos<-cbind(datos,Grupo=grupos$cluster)
head(NDatos) Edad.10 Antiguedad Espacios.Parqueo Velocidad.Cajas
Ariana 2.5 1 7.6 7.6
Guiselle 2.4 0 4.8 9.0
Francisco 2.8 7 6.8 8.4
Griselda 2.3 0 3.4 7.8
Damaris 4.9 6 7.0 3.2
Johana 3.2 4 5.6 7.8
Distribucion.Productos Atencion.Empleados
Ariana 7.8 9.7
Guiselle 7.2 10.0
Francisco 7.6 8.7
Griselda 9.0 10.0
Damaris 1.2 10.0
Johana 6.8 10.0
Calidad.Instalaciones Ubicacion Limpieza Variedad.Productos
Ariana 5.0 9.0 7.6 5.6
Guiselle 2.0 9.6 6.8 8.4
Francisco 2.7 9.2 6.2 9.0
Griselda 1.0 10.0 4.4 4.0
Damaris 4.0 9.0 1.4 4.8
Johana 3.0 10.0 5.0 4.2
Prestigio.Empresa Calidad.Servicio Grupo
Ariana 7.0 6.6 3
Guiselle 9.8 5.4 3
Francisco 9.6 8.5 1
Griselda 2.8 5.4 3
Damaris 2.6 3.3 2
Johana 4.2 7.2 1
# Establezco el directorio en donde se quiere grabar el archivo
#setwd("~/Google Drive/MDCurso/Datos")
# Se graba el archivo en como un CSV
write.csv(NDatos,"NDatos.csv")setwd("~/Google Drive/MDCurso/Datos")
datos <- read.csv("EjemploClientesCorregidaEdad.csv",header=TRUE, sep=";", dec=",", row.names=1)
str(datos)## 'data.frame': 37 obs. of 12 variables:
## $ Edad.10 : num 2.5 2.4 2.8 2.3 4.9 3.2 2.6 2.3 2.5 2.9 ...
## $ Antiguedad : int 1 0 7 0 6 4 0 4 4 0 ...
## $ Espacios.Parqueo : num 7.6 4.8 6.8 3.4 7 5.6 6.2 5.6 4.6 5.4 ...
## $ Velocidad.Cajas : num 7.6 9 8.4 7.8 3.2 7.8 8 6.8 8 6.4 ...
## $ Distribucion.Productos: num 7.8 7.2 7.6 9 1.2 6.8 6.6 6.2 3.8 8.8 ...
## $ Atencion.Empleados : num 9.7 10 8.7 10 10 10 9.3 9.7 10 9.7 ...
## $ Calidad.Instalaciones : num 5 2 2.7 1 4 3 3.3 4 1.7 6.7 ...
## $ Ubicacion : num 9 9.6 9.2 10 9 10 8.6 6.8 9.8 10 ...
## $ Limpieza : num 7.6 6.8 6.2 4.4 1.4 5 7.8 6.8 5 5.6 ...
## $ Variedad.Productos : num 5.6 8.4 9 4 4.8 4.2 6.4 7.4 4.4 6.2 ...
## $ Prestigio.Empresa : num 7 9.8 9.6 2.8 2.6 4.2 9.6 5.6 6 8.4 ...
## $ Calidad.Servicio : num 6.6 5.4 8.5 5.4 3.3 7.2 6.5 4.5 7.6 6.5 ...
grupos<-kmeans(datos,centers=3,iter.max=200,nstart=200)
grupos$cluster # Cluster al que pertenece cada fila## Ariana Guiselle Francisco Griselda Damaris Johana
## 3 3 1 2 2 2
## Bernal Freddy Estafania Laura Arnoldo Beatriz
## 3 1 2 3 1 1
## Rebeca Sofia Ingrid Rocio Karen Luis
## 3 1 1 2 1 3
## Pedro Lorena Elena Julian Natalie Shirley
## 3 2 1 2 2 2
## Andres Alejandro Grace Nuria Flor Roberto
## 2 1 2 3 3 1
## Victor Arturo Maritza Diana Juan Guillermo
## 3 1 3 2 2 3
## Silvia
## 3
grupos$centers # Centros de los clústeres## Edad.10 Antiguedad Espacios.Parqueo Velocidad.Cajas
## 1 3.118182 6.181818 6.527273 8.018182
## 2 3.169231 2.769231 5.400000 7.461538
## 3 2.392308 1.307692 6.492308 8.230769
## Distribucion.Productos Atencion.Empleados Calidad.Instalaciones
## 1 6.927273 9.554545 3.845455
## 2 4.476923 9.692308 2.084615
## 3 8.261538 9.669231 4.284615
## Ubicacion Limpieza Variedad.Productos Prestigio.Empresa Calidad.Servicio
## 1 8.509091 6.818182 7.200000 7.618182 4.781818
## 2 9.200000 5.646154 5.723077 5.261538 5.069231
## 3 9.215385 7.615385 7.676923 8.430769 5.423077
# Gráfico Tipo Araña
centros<-grupos$centers
rownames(centros)<-c("Cluster 1","Cluster 2","Cluster 3")
centros<-as.data.frame(centros)
maximos<-apply(centros,2,max)
minimos<-apply(centros,2,min)
centros<-rbind(minimos,centros)
centros<-rbind(maximos,centros)
centros## Edad.10 Antiguedad Espacios.Parqueo Velocidad.Cajas
## 1 3.169231 6.181818 6.527273 8.230769
## 11 2.392308 1.307692 5.400000 7.461538
## Cluster 1 3.118182 6.181818 6.527273 8.018182
## Cluster 2 3.169231 2.769231 5.400000 7.461538
## Cluster 3 2.392308 1.307692 6.492308 8.230769
## Distribucion.Productos Atencion.Empleados Calidad.Instalaciones
## 1 8.261538 9.692308 4.284615
## 11 4.476923 9.554545 2.084615
## Cluster 1 6.927273 9.554545 3.845455
## Cluster 2 4.476923 9.692308 2.084615
## Cluster 3 8.261538 9.669231 4.284615
## Ubicacion Limpieza Variedad.Productos Prestigio.Empresa
## 1 9.215385 7.615385 7.676923 8.430769
## 11 8.509091 5.646154 5.723077 5.261538
## Cluster 1 8.509091 6.818182 7.200000 7.618182
## Cluster 2 9.200000 5.646154 5.723077 5.261538
## Cluster 3 9.215385 7.615385 7.676923 8.430769
## Calidad.Servicio
## 1 5.423077
## 11 4.781818
## Cluster 1 4.781818
## Cluster 2 5.069231
## Cluster 3 5.423077
color <- c("#EEE6AB","#C5BC8E","#5E412F")
radarchart(as.data.frame(centros),maxmin=TRUE,axistype=4,axislabcol="slategray4",
centerzero=FALSE,seg=8, cglcol="gray67",
pcol=color,plty=1,plwd=5,title="Comparación de clústeres")
legenda <-legend(1.5,1, legend=c("Cluster 1","Cluster 2","Cluster 3"),
seg.len=-1.4,title="Clústeres",pch=21,bty="n" ,lwd=3, y.intersp=1,
horiz=FALSE,col=color)setwd("~/Google Drive/MDCurso/Datos")
datos <- read.table('SpamData.csv', header=TRUE, sep=';',dec='.')
datos<-datos[,c(2,4,5,6,7,9,10,11)]
t1<-system.time(grupos<-kmeans(datos,centers=3,iter.max=3,nstart=1))
t1 user system elapsed
0.004 0.001 0.004
par(mfrow=c(1,2))
centros<-grupos$centers
rownames(centros)<-c("Cluster 1","Cluster 2","Cluster 3")
centros<-as.data.frame(centros)
maximos<-apply(centros,2,max)
minimos<-apply(centros,2,min)
centros<-rbind(minimos,centros)
centros<-rbind(maximos,centros)
color <- c("#FCEBB6","#78C0A8","#5E412F")
radarchart(as.data.frame(centros),maxmin=TRUE,axistype=4,axislabcol="slategray4",
centerzero=FALSE,seg=8, cglcol="gray67",
pcol=color,plty=1,plwd=5,title="Comparación de clústeres")
legenda <-legend(1.5,1, legend=c("Cluster 1","Cluster 2","Cluster 3"),
seg.len=-1.4,title="Clústeres",pch=21,bty="n" ,lwd=3, y.intersp=1,
horiz=FALSE,col=color)
t2<-system.time(grupos<-kmeans(datos,centers=3,iter.max=200,nstart = 200))
t2 user system elapsed
0.388 0.017 0.406
centros<-grupos$centers
rownames(centros)<-c("Cluster 1","Cluster 2","Cluster 3")
centros<-as.data.frame(centros)
maximos<-apply(centros,2,max)
minimos<-apply(centros,2,min)
centros<-rbind(minimos,centros)
centros<-rbind(maximos,centros)
radarchart(as.data.frame(centros),maxmin=TRUE,axistype=4,axislabcol="slategray4",
centerzero=FALSE,seg=8, cglcol="gray67",
pcol=color,plty=1,plwd=5,title="Comparación de clústeres")
legenda <-legend(1.5,1, legend=c("Cluster 1","Cluster 2","Cluster 3"),
seg.len=-1.4,title="Clústeres",pch=21,bty="n" ,lwd=3, y.intersp=1,
horiz=FALSE,col=color)generate <- function(n=50, extradim=0, sigma=1, mu=7) {
data1 <- matrix(rnorm(n*2), n, 2) * sigma
data1[,1] <- data1[,1] + centers[1,1] * mu
data1[,2] <- data1[,2] + centers[1,2] * mu
data2 <- matrix(rnorm(n*2), n, 2) * sigma
data2[,1] <- data2[,1] + centers[2,1] * mu
data2[,2] <- data2[,2] + centers[2,2] * mu
data3 <- matrix(rnorm(n*2), n, 2) * sigma
data3[,1] <- data3[,1] + centers[3,1] * mu
data3[,2] <- data3[,2] + centers[3,2] * mu
data <- rbind(data1,data2,data3)
if (extradim > 0) {
noise <- matrix(rnorm(3*n*extradim)*sigma, 3*n, extradim)
data <- cbind(data, noise)
}
return(data)
}centers <- matrix(c(0,3,1,3,0,4), 3, 2, byrow=T)
centers [,1] [,2]
[1,] 0 3
[2,] 1 3
[3,] 0 4
Data1 <- generate(extradim=0)
head(Data1) [,1] [,2]
[1,] 0.6413689 21.26642
[2,] 0.7323902 21.49162
[3,] 0.1251402 20.46348
[4,] 1.5994146 22.12053
[5,] -0.4862423 18.69831
[6,] 0.1052660 21.79795
dim(Data1)[1] 150 2
groups<-kmeans(Data1,centers=3,iter.max=100,nstart=50)plot(Data1)
points(groups$centers,pch=19,col="blue",cex=2)
points(Data1,col=groups$cluster,pch=19)mydata <- data.frame( (Data1), groups$cluster)
clusplot(mydata, groups$cluster, color=TRUE, shade=TRUE,
labels=2, lines=0)library(scatterplot3d)
Data2 <- generate(extradim=1)
g3d<-scatterplot3d(Data2)
groups<-kmeans(Data2, 3)
g3d$points3d(groups$centers,pch=19,col="blue",cex=2)
g3d$points3d(Data2,col=groups$cluster,pch=19)Datos <- generate(extradim=0)
InerciaIC<-rep(0,30)
for(k in 1:30) {
grupos<-kmeans(Datos,centers=k,nstart=25)
InerciaIC[k]<-grupos$tot.withinss
}
plot(InerciaIC,col="blue",type="b")mis.datos <- scale(Datos)
fviz_nbclust(mis.datos, kmeans, method = "gap_stat")fviz_nbclust(mis.datos, kmeans, method = "wss")fviz_nbclust(mis.datos, kmeans, method = "silhouette")setwd("~/Google Drive/MDCurso/Datos")
Datos <- read.csv("EjemploClientesCorregidaEdad.csv",header=TRUE, sep=";", dec=",", row.names=1)InerciaIC<-rep(0,30)
for(k in 1:30) {
grupos<-kmeans(Datos,centers=k,nstart=25)
InerciaIC[k]<-grupos$tot.withinss
}
plot(InerciaIC,col="blue",type="b")mis.datos <- scale(Datos)
fviz_nbclust(mis.datos, kmeans, method = "gap_stat")fviz_nbclust(mis.datos, kmeans, method = "wss")fviz_nbclust(mis.datos, kmeans, method = "silhouette")setwd("~/Google Drive/MDCurso/Datos")
datos <- read.csv("EjemploClientesCorregidaEdad.csv",header=TRUE, sep=";", dec=",", row.names=1)
grupos<-kmeans(datos,centers=2,iter.max=200,nstart=200)
# Gráfico Tipo Araña
centros<-grupos$centers
rownames(centros)<-c("Cluster 1","Cluster 2")
centros<-as.data.frame(centros)
maximos<-apply(centros,2,max)
minimos<-apply(centros,2,min)
centros<-rbind(minimos,centros)
centros<-rbind(maximos,centros)
color <- c("red","blue")
radarchart(as.data.frame(centros),maxmin=TRUE,axistype=4,axislabcol="slategray4",
centerzero=FALSE,seg=8, cglcol="gray67",
pcol=color,plty=1,plwd=5,title="Comparación de clústeres")
legenda <-legend(1.5,1, legend=c("Cluster 1","Cluster 2"),
seg.len=-1.4,title="Clústeres",pch=21,bty="n" ,lwd=3, y.intersp=1,
horiz=FALSE,col=color)